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Abstract. We present an elementary method to obtain Green's functions in non-perturbative quantum field 
theory in Minkowski space from calculated Green's functions in Euclidean space. Since in non-perturbative 
field theory the analytical structure of amplitudes is many times unknown, especially in the presence of 
confined fields, dispersive representations suffer from systematic uncertainties. Therefore we suggest to 
use the Cauchy-Riemann equations, that perform the analytical continuation without assuming global 
information on the function in the entire complex plane, only in the region through which the equations 
are solved. We use as example the quark propagator in Landau gauge Quantum Chromodynamics, that is 
known from lattice and Dyson-Schwinger studies in Euclidean space. The drawback of the method is the 
instability of the Cauchy-Riemann equations to high-frequency noise, that makes difficult to achieve good 
accuracy. We also point out a few curiosities related to the Wick rotation. 
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1 The Wick rotation 
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Central to non-perturbative quantum field theory is the 
computation of Green's functions, the vacuum expecta- 
tion value of quantum operators. These and the related 
Scattering Matrix elements are most often computed in 
Euclidean space, defined by the transformation 



t 



-He fco — * ik-o 



(1) 



This coordinate transformation is known as "Wick rota- 
tion" (see for example pQ for a short account). There are 
many advantages in solving the field equations in terms of 
the rotated variables to obtain so called Schwinger func- 
tions, and we list some below in section [2] 

Once the wanted functions have been computed in Eu- 
clidean (momentum) space Re — (ifco,k) one would wish 
to recover the original Minkowski space Green's functions 
by inverting the Wick rotation. This is possible in pertur- 
bation theory at low orders where one has explicit 
expressions for the functions, and their analytical struc- 
ture (poles, cuts, essential singularities) is at hand, so that 
one can employ Cauchy's theorem and collect dispersive 
cut integrals or pole residues if need be, and obtain the 
Minkowski space Green's function by analytical continua- 
tion. 

For example, to obtain the electron propagator in mo- 
mentum space (Fourier transform of the probability am- 
plitude for the electron to reach point x if it was originally 



S(p) 



d A xe- ixp S(x,0) = iZ{ P 2 )- 



M(p 2 ) 



■ p z _ M 2 (p 2 ) +ie 

one would perform the Wick rotation po — > ipoE and ob- 
tain a function S(j>e) as a perturbation of its free-field 
values (M(p 2 ) = m, Z(p 2 ) = 1). With the function ex- 
plicitly known, one just extends it into the complex plane 
and simply substitute its argument pqe by — ipo- If the 
function is not totally known but its analytical structure 
is, one employs Cauchy's theorem as mentioned. 

However, in non-perturbative quantum field theory one 
is seldom in this desirable situation. More often than not, 
the function has been calculated with the help of a com- 
puter, be it by solving the Dyson-Schwinger equation [4] 
with some carefully designed truncation, or by trying a 
Montecarlo evaluation averaging over a small number of 
configurations on a lattice [5]. 

The outcome is that the function is then known for 
Euclidean momenta, typically k% > 0, and an extension 
into the complex plane becomes necessary to reach the 



negative axis k E < 0, 



-k . For non-perturbative 



functions one sometimes ignores the precise analytic struc- 
ture in the complex plane. This situation is worsened in 
theories where the field quanta do not appear in asymp- 
totic states, except in very specific combinations, such as 
is presumably the case for Quantum Chromodynamics. 

Attempts have of course been made to solve the prob- 
lem. An obvious approach is to perform a theory-motivated 
fit to the computer data, building-in well educated guesses 
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Fig. 1. The analytical continuation to negative Euclidean 
squared momentum (positive quark virtuality in Minkowski 
space) of the Landau gauge quark mass function. The yel- 
low band is obtained with input from an Euclidean Dyson- 
Schwinger calculation, where the Cauchy-Riemann equations 
have been solved with the 8 method in different grids. The 
Lattice data is likewise analytically continued with the Cauchy- 
Riemann equations. At low p 2 the error is dominated by the 
original lattice statistical error, eventually the build-up of nu- 
merical errors in the Cauchy-Riemann equations dominate. 
The trend of Re(M) is clearly decreasing with p 2 , and therefore 
a pole in the quark propagator is expected. Since Im(M), pre- 
sented later on, is small, this pole is on or close to the real axis, 
with M = 305(25) MeV. The lattice data that is analytically 
continued is from ref. [11] , 



on what the analytical structure of the continued function 
must be [g|Tf]. 

Another possibility is to write a spectral representation 
in terms of a Stieltjes transform yielding a spectral density 

P, 

p(m 2 ) 



S(P 2 ) = J dp~ 2 



(3) 



(we have eliminated spin structure) and then solving the 
Dyson-Schwinger equation directly "in Minkowski space" 
for p, finally inverting the Stieltjes transform to recover 
the propagator or Bethe-Salpeter wanted function [8], [9]. 
It is of course clear that in writing the spectral representa- 
tion, one is already assuming a given cut structure for the 
function. It is known that, for local quantum field theories 
where the quanta appear in the final state, the propaga- 
tor in coordinate space in the upper-half x complex plane 
is analytic [3J, but for confined quanta many questions 
remain. 

In this paper we enrich the toolbox by putting forward 
a very simple method that does not require the function's 
analytical structure to be known on the entire complex 
plane. We observe that the analyticity of a given function 
cannot only be formulated globally, through satisfaction of 
Cauchy's theorem, but also locally, through satisfaction of 
the Cauchy-Riemann equations. By integrating this sim- 



ple first order differential system with initial condition the 
function in a given region computed by other means (DSE, 
lattice, exact renormalization group equations [10j . 
etc.) one can achieve two goals. First, the numerically ob- 
tained solution can break down at a given point or line in 
the plane, indicating perhaps a pole or other singularity. 
Second, if the region where the system is integrated avoids 
such singularities, one can obtain the analytical extension 
(within errors and non-uniqueness) to another region in 
the complex plane. 

We use as example the behavior of the quark mass 
function in Landau gauge QCD, M(p 2 ), that we analyti- 
cally continue from positive Euclidean virtuality (p 2 > 0) 
into Minkowski space with (p 2 < 0). This is plotted in 
figure [T] It can be seen that, within the statistical errors 
inherited from lattice data and the systematic numerical 
errors intrinsic to our procedure, the mass function de- 
creases with increasing Minkowski p 2 . Our result keeps 
open the possibility of a pole of the quark propagator for 
real p 2 (whose absence has been at times thought as a 
possible sign of confinement). 

The rest of the paper consists of four sections. In sec- 
tion [2] we make a few comments, some common place 
but others quoted less often, about the advantages of ini- 
tially working in Euclidean space. In section [3] we present 
the Cauchy-Riemann method with one practical case, the 
quark propagator. A few theoretical comments about er- 
rors involved in the process and the generalization to more 
dimensions are left for section 01 Our discussion is sum- 
marized in section [S] 



2 Working in Euclidean space 

In lattice formulations of Quantum Field Theory, the field 
configurations over which the path integral is evaluated 

are randomly generated according to a distribution e~ f dtC , 
that is the Wick-rotation of the actual quantum weight 

for the path integral, e ■/" dtC . Thereafter computed lat- 
tice Green's functions are valid in Euclidean space. Even 
if working in Minkowski space, a popular way of "min- 
imally" regularizing in the path integral formalism is to 
rotate the time integration into the complex plane 



d x = lim 



dt 



Beyond the convergence of the path integral and the 
weighting configurations inside a compact set of function 
space, there are several more advantages. 

One is that the Dyson-Schwinger equations, whose so- 
lutions are often used to interpret lattice data, are ex- 
tremely difficult to solve on a computer in Minkowski 
space. Indeed, a typical DSE is that for the mass func- 
tion of a fermion in the presence of a scalar field, with 
Yukawa coupling, in rainbow approximation 



M(p) 



d 4 k 



M(k) 



{k 2 - M(k) 2 + ie)((k - p) 2 -m 2 + ie) 



(4) 
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If one solves the equation iteratively by guessing Mq , one 
notices that the position of the fermion pole in the denom- 
inator is not known after the first iteration, and it needs 
to be determined numerically (an attempt at carrying on 
this program exists [T!?]). 

However the common use is to Wick-rotate the k inte- 
gration variable to Euclidean space. If p is likewise rotated, 
the resulting equation is easier to program as the denom- 
inator poles are on the left k 2 plane, out of the numerical 
integration region in the radial k 2 variable. 



M(p) = cjd i k 



M(k) 



{k 2 + M{k) 2 + ie)((k~p) 2 + m 2 + ie) 



(5) 

(c represents constants irrelevant to the discussion). 

Note also that when working in Euclidean space, a dis- 
crete subgroup of the Euclidean rotation group is retained. 
For example, for a simple cubic lattice with (x = ai x , y = 
ai y , z — ai z , t = ait), invariance under rotations by tt/2 is 
explicit. This can be exploited to study the quantum rep- 
resentations of the discrete group, then trying to match 
the resulting states to a representation of the full contin- 
uous group. 

However, in Minkowski space there is no finite lattice 
that retains invariance under a non-trivial subgroup of 
the Lorentz group. If a grid is invariant under discrete 
Lorentz transformations of parameter a, then it has in- 
finitely many points. We discuss Lorentz-invariant dis- 
cretizations of Minkowski space in the appendix. 

One further motivation is the non-compactness of the 
equal k 2 hypersurfaces. While in Euclidean space the con- 
dition k 2 = A 2 determines a hypersphere's surface, so that 



d 4 k E f(k 2 E ) = 2tt 2 



k%dk E f(k 



(6) 



can be factorized into a radial integral and a finite 27r 2 
solid hyperangle (the hyperarea of a unit-radius hyper- 
sphere's surface) this is not possible in Minkowski space, 
where the corresponding unit-hyperboloid k 2 = 1 has infi- 
nite hypersurface. Therefore, integrals of Lorentz-invariant 
functions are by necessity divergent even after regulation 
of large virtualities, and are only defined by analytical 
continuation from Euclidean space. 

In perturbation theory, a much used method is to per- 
form the ko integrals first, usually with pole analysis, and 
later impose a cutoff on spacelike momentum k. How- 
ever this cutoff is frame-dependent, and there is no direct 
method that manifestly preserves Lorentz invariance. 

From all these arguments, it is hard to conceive progress 
in non-perturbative quantum field-theory in Minkowski 
space without progress in complex-plane analytical contin- 
uation for the relevant functions. This paper is a modest 
contribution in this direction, with the interest in keeping 
the discussion alive. 



3 Numerical solution of the Cauchy-Riemann 
equations 

3.1 Cauchy-Riemann equations in polar coordinates 

If u and v are respectively the real and imaginary parts 
of a complex function of one complex variable p 2 = re , 
the Cauchy-Riemann equations in polar coordinates read 



dv du 
~d() =r ~dr~ 
du dv 
06 ~ ~ T d0 ' 



(7) 
(8) 



Given the initial conditions it(r, 0) = Uo(r), v(r, 0) = 
vo(r) on a segment of the realp 2 axis, corresponding to 9 = 
0, one can then evolve the system towards increasing and 
decreasing 9 (like the opening of a fan). For very smooth 
data sets one can typically reach 90-120 degrees on each 
side of the fan before the instabilities wipe the solution to 
infinity. The (Cauchy-Euler) explicit discretization with 
centered r-derivative on a grid (r 3 , 9{) is simply 



v(r 3 ,9 i+1 ) = v(r 3 ,9i) + r 3 (Q i+1 - $i) 



u(r 3+1 ,9i) - u(rj-i,6i) 



r j+i - rj-i 



(9) 



(10) 



where, to solve over an arch taken anticlockwise, (#;+i — 
6j) > 0. At the end-points of the grid one cannot use 
centered derivative, so left (right derivative) is necessary, 
{v(r2, 9i) — v(r\, 9i))/{ri — r\), etc. The situation is repre- 
sented in figure [2] 




Fig. 2. The Cauchy-Riemann equations in polar coordinates 
allow to explore a fan-shaped region of the complex plane 
where a function is analytic. 



The use of this method is to provide a cross-check of 
solutions of the Dyson- Schwinger equations in the com- 
plex plane. These are needed to solve the Bethe-Salpeter 
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equations for mesons, since the (external) meson momen- 
tum is of course in Minkowski space (real), and the inter- 
nal quark momentum is Wick-rotated to Euclidean space 
(imaginary) so that one ends solving the DSE inside a 
parabolla in the complex plane symmetric respect to the 
real momentum axis. The Cauchy-Riemann equations are 
currently no match in precision to directly solving the DSE 
in the complex plane where this is feasible, but they can 
provide a cross check that is very simple to programme 
(compare the trivial linear system above with the complex, 
non-linear, bidimensional DSE when the angular kernel or 
vertex are non-trivial). 

The Cauchy-Riemann equations however are a state- 
ment of analyticity, and the solution is a numerical repre- 
sentation of the closest analytical function that contains 
the initial data. This means that if the "true" function 
has a pole or a cut, the Cauchy-Riemann iteration will 
fail to see it and simply separate from that function, and 
is likely to diverge soon from accruing instabilities. This 
is illustrated in figure O 



3.2 Cauchy-Riemann equations on a strip 

The advantage of a local formulation of analyticity em- 
ploying the Cauchy-Riemann equations is lost if one needs 
to swipe the entire complex plane. Therefore it is prof- 
itable to solve them in Cartesian coordinates, first away 
from the x axis along y 



du dv 

dy dx 

dv du 

dy dx 



(11) 
(12) 
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and then leftwards along x. The method fails if the com- 
plex plane is completely cut from — oo to oo along the y 
axis. In any other situation (the standard half-plane cut 
of a power-law or logarithm, or a finite number of poles 
or essential singularities) , one can find a path between the 
right and the left x axis and solve the Cauchy-Riemann 
equations along them. 

We now improve upon the discretization of the differ- 
ential equations and employ an implicit 9 method. This is 
convenient since the Cauchy-Riemann are quite unstable 
(as pointed out below in subsection 14. 2p . The d y u equa- 
tion for a point not on the edge of the grid becomes (with 
the superindex labelling y, the subindex x) 



if 



!/.' 



1 



(13) 



In the advance along y one groups the ui and Vi for 
fixed yj in a vector u J = (u{, v\, u^, v%, u J N , v J N ) and the 
9 method's discretization can be written down as a linear 
problem 

Au J+1 = Bu 3 . (14) 



for z — 



Fig. 3. Real and imaginary part of the function 
re 10 with fixed r = 0.99 and varying 9. The dashed line is the 
computer solution of the Cauchy-Riemann equations in degree 
steps, with the initial condition given on the positive real half- 
axis. The function has a pole at (1 + i)/\/2 that the Cauchy- 
Riemann equations cannot isolate since they entail analyticity. 
They however diverge due to accruing instability in the region 
where the function has larger derivative. 



In the simplest case of equal x subintervals, one can 
define r = 9Ay/(2Ax), then the matrix A becomes 



f 1 


-2r 





2r 





2r 


1 


-2r 










— r 


1 








r 








1 


— r 





2r 



-2r 




1 2r 

V ... 2r -2r 1 / 

where the third and fourth line are the repeated unit, ex- 
cept the non-vanishing elements are shifted to the right 
(the l's always mark the diagonal, and the matrix is band- 
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Fig. 4. The Cauchy-Riemann equations in Cartesian coordi- 
nates can be used to explore a strip, first upwards along the 
j/axis, then leftwards along the x axis. In this use one only 
needs analyticity of the function on a strip above (or below) 
the axis to obtain information on the Minkowski side. The 
method fails only if the complex plane is completely cut from 
— oo to oo along the y axis. 



diagonal) . The matrix B can likewise be filled by exchang- 
ing 9 — » (1—9) and changing the sign of all off-diagonal 
matrix elements of A (the diagonal of B likewise contains 
1's). 

The advance in y proceeds by solving the linear system 
to obtain (u, v) at yj+i from their values at yj. We use 
standard LU factorization, although since the matrix A is 
band-diagonal, Crout's algorithm can speed things some. 
For one complex variable and the small number of points 
we use this is irrelevant in computer time. 

However, since the advance to the left will have as 
initial condition the edge (along the y axis) of the first 
computed block, and at the end we will be interested on 
the values of the function on the x axis, that is the edge 
of the second computed block, it pays off to improve the 
computation of the derivative at the edge of the block. 

The left-derivative we have displayed explicitly, 



du(x — 0, y) u J 2 



o(h) 



(16) 



dx x-i — x\ 

can be interpreted as a centered derivative at the mid- 
point (x2 + x\)/2. Considering also the centered deriva- 
tive at X2 and extrapolating linearly to x\, one obtains an 
improved 



du(x = 0, y) 
dx 



x 2 — X\ x-i - x\ 



o(h 2 



(17) 



brief eigenvalue analysis is presented below that explains 
why, in subsection 14.21 

To show a test of the method, we employ the simple 
function 

W*) = YTz^l (18) 

that has two poles above and below the x axis at ±2z. 
Because of the fast build-up of numerical errors, we need 
Ay < Ax so the strip is always shorter in the direction of 
the advance of the integration (for an N x N problem). 
In figs. [5] and [6] we show the real and imaginary parts cal- 
culated with the Cauchy-Riemann equations with initial 
condition on the positive real half-axis, plotted along the 
imaginary and the negative real-axis respectively. 



Method check 



and likewise at the last N point of the grid, and this 
slightly complicates the first two and last two rows of A. 

Once the advance upwards in the y direction has reached 
N, one starts an advance to the left as in figure|4]and sim- 
ilar considerations apply. 

The 9 parameter that advances or delays the deriva- 
tive perpendicular to the integration direction is empiri- 
cally fixed for now. Several problems are somewhat inde- 
pendent of 9, others have a broad minimum of instability 
around -1. We find that 9 — —0.5 is as good as any. A 
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Fig. 5. The function /test(«) = 1+^/4 analytically continued 
from (x > 0, y = 0) to (x — 0, y > 0) with the Cauchy-Riemann 
equations. The imaginary part is exactly zero on this imaginary 
axis, the function has a pole at y = 2. 



The method starts breaking down for x near zero and 
y > 1 as the pole is approached. But one can see how 
the part of the strip that goes well below the pole passes 
cleanly and allows to reproduce the function on the left 
axis, given as initial condition the N exact values on the 
right axis. All in all, the method provides a reasonable 
representation of the function from the solution of the 
Cauchy-Riemann equations. Note the imaginary part, ex- 
actly zero on the left axis, is calculated to be of order 1% 
with a forty-point grid. This should be considered the er- 
ror of the method, and is far less than the statistical errors 
in the lattice data that we will shortly employ. 

Should the analytic structure of the function become 
available, one could devise an arbitrary path in the com- 
plex plane from the region where the function is known to 
the region where it is wanted by analytical continuation. 
The initial value problem can then be formulated with the 
advance direction along the tangent vector to the path, r. 
This vector changes in principle orientation, so one would 
need to use the "Cartesian-like" formulation to advance 
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Fig. 6. The function /test (2) = 1+ ^ / / 4 analytically continued 
from [x > 0, y = 0) to (a; < 0, y — 0) with the Cauchy-Riemann 
equations. 
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Fig. 7. As in figure [5] but for the function cos z. 



and the "polar-like" formulation to rotate the direction of 
advance, in alternate steps. 



3.3 Wick rotation of the quark mass function 

In the introduction we have advanced, in figure[l] our main 
application, the analytical continuation of the quark mass 
function M(p 2 ) in Landau gauge, that we now carefully 
analyze and justify. 

Because of Lorentz invariance applied to a spin 1/2 
fcrmion, the quark propagator can be written in full gen- 
erality as 

iZ(p 2 ) £> iZ(p 2 )M(p 2 ) 
{P> p 2 -M 2 (p 2 ) p 2 -M 2 (p 2 ) ' 1 ' 



Fig. 8. As in figure [6] but for the function cos z. 

As already discussed, the pole in the denominator is a nui- 
sance usually disposed of by performing the Wick rotation 
Po — > ip%- Concentrating on the denominator then 



1 



-1 



Pi 



M 2 (p 2 



vl 



p 2 + M 2 (-p 2 ~p 2 ) 



the pole is absent for real M. One usually eschews a sign, 
the function M 2 (— p 2 , — p 2 ) is customarily called M(p 2 E ) 
where p 2 E > 0, and we will keep this notation. Hence to 
retrieve the mass function that actually appears in the 
Minkowski space propagator for positive virtuality p 2 , we 
need to identify it with the analytical continuation of the 
lattice (or DSE) M(p 2 E ) to negative p 2 . Note that although 
the analytical continuation is nominally made inp 0l since 
Pq — p 2 is a polynomial, it is an analytical function of 
Po . We conclude that wherever M is analytical in po it is 
also analytical in p\ — p 2 (conversely, there needs to be a 
branch cut in the po plane that is absent in terms of the 
Lorentz invariant variable). 

The analytical properties of M(p 2 ) are not well known, 
especially for confined quanta such as quarks. However we 
note that a pole in this function would imply a zero in the 
quark propagator, and this, assumed continuous, can be 
ruled out in the entire region of the complex plane that 
we sample, from the lattice data on the Euclidean real p 2 
axis. 

We take the lattice data from [TT] . This has been pro- 
vided to us from 28 3 x 96 lattices with MILC configura- 
tions. We have used a -1 = 2.29 GeV to set the scale from 
the internodal spacing. 

Due to our sensitivity to large frequency noise, we only 
use a subset of the lattice data, taking one of every few 
points in the interval p G (0, 8) GeV. The trimming has 
been performed so that the resulting mass monotonously 
decrease towards higher momenta (asymptotic freedom), 
to avoid distortions of analyticity and large errors through 
rapidly varying derivative. We further square the abscissa 
p — > p 2 , since the latter is the variable in which we per- 
form the analytical continuation. The original lattice data, 
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Fig. 9. A part of the lattice data for the quark mass func- 
tion in Landau-gauge QCD. We have trimmed the data to en- 
sure the monotonous decrease in the function and reduce the 
high-frequency noise, which grows fast in the Cauchy-Riemann 
equations. We slightly increased the error bands to cover the 
omitted points (trading our systematic error in trimming into 
a statistical error). The lattice data from [11] is normalized to 
a current quark mass m u = 60 MeV. 



Fig. 10. Input calculation from the Dyson-Schwinger formal- 
ism. The current quark mass is 2 MeV at 13 GeV. The function 
is smoother than lattice data, reducing high-frequency noise, 
in exchange its systematic errors are more difficult to control. 
The original function reported in [16] decreases at very small 
momentum, we have been conservative and avoided this by a 
small variation of the vertex dressing functions. 



already in physical GeV units, is given in figure [HI We 
also show the actual input set to our code after these ma- 
nipulations have been performed, that is faithful to the 
original data and error bands, but amenable to analyti- 
cal continuation. The data is renormalized in the MOM 
scheme, and as can be seen the mass at a scale of 8 GeV 
is 60 MeV with the chosen scale a -1 . We further take as 
input a Dyson-Schwinger calculation from [16j . This data 
has as an advantage that there are no statistical errors and 
the function is very smooth. In exchange, there are sys- 
tematic errors (coming from the precise way in which the 
quark-gluon vertex is treated in that reference), that are 
unknown and only controllable in the propagator in com- 
parison with lattice data or renormalization group equa- 
tions. We plot the resulting set in figure ITTJl 

Finally, we perform the analytical continuation on a 
strip in the complex plane above the axis, that presum- 
ably avoids non-analyticities in M(p 2 ) (else a continuation 
under the axis is possible) and obtain the real part of M 
advanced in the introduction in figure [T] 

We also plot in figure [11] the imaginary part of the 
same M function, that as can be seen is compatible with 
zero within the error bands. 

From the graphs one can conclude that, just as for 
the tree-level propagator in perturbation theory, there is 
a crossing of M and p for (negative Euclidean), positive 
Minkowski p 2 . This means that the actual quark propa- 
gator does have a pole at or very near the real axis. It 
has been quoted [6] at (300 — 500) M eV from the Dyson- 
Schwinger equations alone. From the analysis of the lattice 
data set at hand, we conclude M(M) = 305(25) MeV, in 
agreement with that estimate. Of course, it would be in- 



Quark mass in Minkowski space 

Landau gauge 
0.2, , 1 , , ^- 



0.15 - 
0.1- 
> 0.05- 




Fig. 11. The Cauchy-Riemann method leads to an imaginary 
part of the mass-function on the left half-axis in Euclidean 
space that is well-compatible with 0. This is in agreement with 
the outcome of a Taylor expansion around the origin (however 
there is no telling the convergence radius of such a series, so 
we deem the Cauchy-Riemann method superior). 



teresting to compare with other lattice data sets, and in 
particular use different current quark masses, so the error 
band is definitely larger. 

We do not find support for the attending conjecture of 
two conjugate poles with a sizeable imaginary part. 
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4 Some theoretical issues 

4.1 Uniqueness 

Here we study to what extent the solution of the Cauchy- 
Riemann equations for the quark mass function is unique, 
given the initial conditions as the lattice computation on 
the positive p 2 E half-axis. By standard complex analysis, 
the uniqueness of an analytic continuation of a function is 
guaranteed if the function is initially known on an open 
subset of C. 

The positive real half-axis is open in R, but not in C. 
However it is easy to show that the analytic continuation is 
unique. Imagine that u and v are known for y = and x > 
0. Then, all partial derivatives ^-^ and are known. In 
particular for the quark mass function, v(y = 0) = (and 
all x-derivatives also vanish), and u(x, y = 0) = M(x), the 
real mass. 

Assume that the extension of M to the complex plane 
was not unique. Then in addition to f — (u, v) there 
would be another function, f + g, that would satisfy the 
Cauchy-Riemann equations with the same initial condi- 
tions. Since the sum of two analytic functions is analytic, 
g itself should be analytic. This means that its components 
g x , g y would also satisfy the Cauchy-Riemann equations 



dg x _ 

dy 
dgy 
dy 



dx 
gx_ 
dx 



(20) 
(21) 



with initial condition g{y — 0) = exactly, with all deriva- 
tives - 9 S y r ^ — also vanishing on the real axis. Auto- 
ax™ ° 

matically, employing cq. (|20|) . and subsequently deriving 
it, all y derivatives also vanish. Therefore g is exactly zero 
in the domain of analyticity, and / unique. 

Of course, in practice / is only known at a discrete and 



finite set of points Zi 



, 0) i = I..N. An analytic func- 



tion could oscillate between any two of the points and take 
arbitrarily large or small values. Therefore one needs an 
additional hypothesis to claim that the computed function 
is a fair representation of the "actual" function. 

The sufficient hypothesis is monotony of the function 
between any two grid points (note the function might 
be globally non-monotonous by being allowed to change 
derivative sign at the grid points themselves). If the func- 
tion is strictly decreasing between Xi and Xi+i, then the 
maximum and minimum values that it can take between 
them are fi and /j+i, and the function is bound (it being 
analytic, it is also continuous). Then, to arbitrarily shrink 
the error in our knowledge of the initial condition, one 
just needs to arbitrarily shrink the grid spacing, so as to 
further constrain the function in every subinterval. The 
function computed with the discretized Cauchy-Riemann 
equations will be as close to the true function as the sta- 
bility of the system allows, given the bound error in the 
initial conditions. 

For our example, the light quark propagator, there is 
essentially no question that the mass function is monotonously 



decreasing towards larger momenta. This is known at large 
momentum from asymptotic QCD and at low momentum 
from all studies of Dyson-Schwinger equations and lattice 
(where all non-monotonous behavior has way less than la 
significance and can be safely called noise) . The hypothe- 
sis of monotony can be checked (falsified) with lattice data 
by simply decreasing the link size in the grid while at the 
same time reducing the statistical error bar. 

One more caveat can be raised. Imagine adding to 
the "actual" function f(u,v) another analytical function 
g(u, v) such that g is very near zero on the right (Eu- 
clidean) axis and very large on the left hand (Minkowski) 
side. Then, while f + g does not exceed the error bars for 
/ on the initial data, it completely changes the answer on 
output since f+g is very different from / on the right half- 
axis. Of course, the derivative of the function must be very 
large around u = since the function changes from very 
small to sizeable in a small interval. To bind this deriva- 
tive from above and exclude this unpleasant possibility 
one needs to demand an additional condition, since exact 
knowledge of all derivatives of the function or knowledge 
of the function in the entire interval is, in a computer grid 
representation, unavailable. 

Now, a fast change of the derivative beyond what is 
visible from the data points implies that the second deriva- 
tive is not well represented by its discrete approximation. 
Here one can demand monotony of the second derivative 
of the function between the last three (few) points of the 
grid on the interval x\, X2, x%. This guarantees that the 
extrapolation of the derivative at just the last grid point 
does not grow arbitrarily, since the second derivative re- 
mains bound. This is now quite a technical condition, and 
maybe not optimal, others being possible. 



4.2 Instability of large frequency noise 

Let us now consider the effect of a perturbation on the 
system of Cauchy-Riemann equations (|20)) . Since the sys- 
tem is linear, it accepts a Fourier analysis. Let us perform 
it on the x variable so that the Fourier components are 



g x {x,y) = A(k,y)e lka 
g y {x,y) = B(k,y)e tkx 

Then the system of equations becomes 

d_ 

dy 



A 
D 





ik 



—ik 




.4 
B 



and therefore 

A 
B 



cosh ky 
i sinh ky 



-i sinh ky 
coshfcy 



A Q 
Bo 



(22) 



(23) 



(24) 



in terms of the initial condition on the x-axis. Obviously, 
if the exact solution is initially perturbed due to computer 
inaccuracies by an amount SAq, with Bq = SBq — for 
simplicity, then at large distances the perturbation on the 
computed solution exponentiates 

5 A , SB cx e kv \SA \ . 
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This could of course be anticipated by remembering 
that the solutions of the Cauchy-Riemann equations are 
two-dimensional harmonic functions, so that the separa- 
ble solutions are sinusoidal functions times exponentials, 
cos kx cosh ky, etc. 

This is especially worrysome when using Montecarlo 
data for the initial condition, since the local (high-fc) noise 
spoils stability very soon. Thus, one needs to apply a 
cooling algorithm or trim the data first to remove short- 
distance fluctuations, justifying our keeping only part of 
the lattice data to ensure monotony. Large-distance, sys- 
tematic shifts of the initial condition are less perturbing. 
In figure [12] we capture the u (real part) noisiest eigenvec- 
tor of the iteration matrix A~ 1 B for a particular 9 = —0.5 
method for fixed grid size, to show its increasing wavenum- 
ber. 



Most unstable eigenvector 
u(x) 



1 1 1 1 1 


■ -N=10 




• ■■■•N=20 




♦ — »N=30 


m 

■ ■ w- 








.** 


~~~~m ■ 




r ' \ 

'• 


A/vVv 


V 



2 4 6 8 10 

x 



Fig. 12. The unstable most eigenvector as a function of the 
number of points transverse to the direction of advance, for the 
iteration matrix in the Cauchy-Riemann 6 — —0.5 discretiza- 
tion. The functions have been vertically shifted for visibility. 



The iteration matrix A _1 B always has eigenvectors 
that are larger than 1 in modulus. We have studied them 
for very simple equispaced rectangular grids, with the di- 
rection of advance along y. For fixed ratio of the incre- 
ment ^f, the largest eigenvalue is quite independent of 
the number of grid points (obvious from the definition of 
the matrix). The dependence with 9 and can be fol- 
lowed from table [T] 



Note from the table that, unlike for the heat equation, the 
9 method is not convergent. We are not able to approach a 
given function with arbitrary accuracy, but only provide 
an estimate. It is apparent that decreasing the advance 
step 2j is n °t a winning strategy, since, for example, fixing 
9 = —0.5, with — 1; 1^1 = 1-6 and to advance the 
same distance with = 1/10, one needs ten steps, but 



Table 1. Largest eigenvalue (one of a pair) for the iteration 
matrix of the 6 method. 



6 

At, 


|A| 


2.5 


1 


9.8 


2.0 


1 


2.6 


1.5 


1 


11 


1.0 


1 


13 


0.5 


1 


2.7 


0.1 


1 


2.0 


-0.1 


1 


1.8 


-0.5 


1 


1.6 


-1.0 


1 


11 


-1.5 


1 


13 


-2.0 


1 


2.7 


-2.5 


1 


7.8 


-0.5 


i 

2 


1.4 


-0.5 


1 

4 


1.2 


-0.5 


1 


1.14 


-0.5 


? 


1.11 


-0.5 


f 

10 


1.08 



1.08 10 = 2.1 > 1.6, meaning that with a smaller step, one 
can advance less far in the progress direction since errors 
amplify faster. (Of course, by decreasing the step one does 
obtain a more reliable representation of the function for 
short advance distances). 

4.3 Generalization to several complex variables 

In principle it would appear straightforward to general- 
ize the Wick rotation to several dimensions. For exam- 
ple, let us consider a vertex function in field theory, say 
the quark and gluon or the electron-photon three-point 
functions. These are characterized by twelve Dirac ten- 
sors multiplied by amplitudes of the three independent 
Lorentz scalar variables, the squared momenta of each of 
the particles, 

Hplplq 2 ) i = 1-12 . 
Given the lattice data in Euclidean space, 

>>i(PiE,P2E><l%) * = 1-12 

one would need to perform the inverse Wick rotation to 
negative p\ in each of the variables, in practice solving the 
Cauchy-Riemann equations variable by variable. Notice if 
the power-law solutions of [16] are correct, then one ex- 
pects a cut at zero virtuality q 2 = in the gluon variable, 
but this can be avoided by appropriately deforming the 
region where one solves the Cauchy-Riemann equation. 

Now let us examine a curiosity that does not come 
about in one dimension. If there is only one variable, 
p 2 = defines a light-cone in Minkowski space, a three- 
dimensional manifold in four-dimensional space. However, 
p 2 E — defines the origin in Euclidean space, just a point 
(this is just another manifestation of the difference be- 
tween the compact rotation group and the unbound Lorentz 
group). The interesting observaion is that, upon Wick- 
rotation, f(p% — 0), the value of a Green's function at one 
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point in Euclidean space, becomes f(p 2 = 0), the value of 
the same Green's function on the entire light-cone. 

But what happens in more dimensions? One may know 
the value of the function \{p 2 E1 — Q,Pe2'Pe2) m Euclidean 
space at the origin for the variable pei- But in Minkowski 
space p\ = does not imply p\ = 0, hence the function 
takes different values for different points of the p\ light- 
cone, \{p\ — 0,^2; (Pi — P2) 2 that do not coincide with the 
value at the origin in Euclidean space. 

This comes about because in a three-point function 
there are two reference four- vectors, p\ and pi, and while 
the Euclidean point with p\ = is at a fixed distance from 
the poin 



4.4 Taylor expansion 

Analytical functions accept Taylor expansions of the type 

n=0 

One could think of performing a polynomial fit to a given 
set of data points (zi,f(zi)) to represent the function 
within the radius of convergence of the series. One would 
simply need to solve the system of Vandermonde for N 
points 

= E(*-W^^) • (26) 

71=0 ^ ' ' 

For example, expanding around the origin, the matrix of 
coefficients for the linear system is the Vandermonde ma- 
trix with rows (1, Zi, z 2 , . . . zf 1 ). Once the system has been 

solved for the derivatives J z °^ , they can be substi- 
tuted in eq. (|25p to obtain the function at an arbitrary 
point. 

One would argue that this is the simplest local im- 
plementation of analyticity, and why should one worry 
about the Cauchy-Riemann equations at all. Of course, 
in practice monotony is difficult to achieve with a finite 
number of polynomials: the approximant will oscillate be- 
tween the tabulated grid points with the lattice data. In 
addition, should there be a cut starting at p 2 = in the 
complex plane, that would not be surprising in view of 
the power-law representations reported in the literature 
for other Green's functions [6], the radius of convergence 
of the Taylor series would be exactly zero. Although such 
difficulties can be circumvented, a practical implementa- 
tion would become as or more difficult than the Cauchy- 
Riemann equations. We have not pursued the matter fur- 
ther. 



5 Summary and conclusions 

We have presented a first analysis of the Cauchy-Riemann 
equations as applied to performing the inverse Wick rota- 
tion from Euclidean to Minkowski momenta. Given that 



their square is the elliptic Laplace equation, and that they 
are not of variational type [17j . we do not have a trivial 
method to solve them to arbitrary accuracy. However they 
provide an estimate of a function by analytical continu- 
ation as an initial value problem, that relies only on the 
existence of a path on the complex plane that is free of 
singularities. While much mathematical literature exists 
[18j and more sophisticated solution methods have been 
reported, the very elementary methods presented here are 
sufficient to gain insight into the quark propagator mass 
function. 

We obtain, after analytical continuation of lattice QCD 
data, qualitatively backed up by a Dyson-Schwinger cal- 
culation, and in agreement with recent studies, a pole of 
the quark propagator in Landau gauge QCD at or very 
near the real axis, with a mass of 305(25) MeV . We find 
no support for the sometimes conjectured pair of complex 
conjugate poles. 
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A 

Lorentz invariant discretizations in Minkowski 
space 

Physical states provide representations of the rotation group 
In lattice gauge theory, the symmetry of the grid reclas- 
sifies the states but at least offers some control over what 
signal may belong to what spin multiplet. This is because 
the lattice is invariant under a subgroup of the rotation 
group, typically tagged by a minimum angle ~ 90 de- 
grees. Likewise, the grid is invariant under a subset of the 
translation group characterized by the grid spacing a. 

However there is no equivalent parameter for finite 
grids in Minkowski space. A grid that is invariant un- 
der discrete Lorentz transformations has infinitely many 
points (and is therefore not tractable on a computer). 
We discuss this shortly as one more advantage of working 
in Euclidean space that can come to play when comput- 
ing pdf 's or reducing Green's functions in non equal-time 
gauges such as the light-front gauge. 

To see it, it is easiest to observe that Lorentz transfor- 
mations map the light-cone to itself, acting as dilatations 
within this manifold. To see it, it is simpler to work in 
1+1 dimension, where the light-cone becomes the pair of 
lines 
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The action of the Lorentz boost becomes 



tx 



7 



7 



to 

x = t Q 



hence 



h = t (7 + fa) , 



(27) 



(28) 



a dilatation of parameter 7(1+/?). There is now an obvious 
way to construct a discretization of the light-cone that is 
invariant under Lorentz transformations of a fixed param- 
eter 7. Simply pick an arbitrary point (to,xo) and obtain 
the sequence (ti,x±), (^a^)--- obtained by successively 
applying eq. (f2"T|) to it. The infinite sequence so obtained 
is such that every point j is the image of another point 
j — 1 under a Lorentz transformation, except for j = 0. 
To generate this one we further need to also include in the 
discretization all points obtained by successively applying 
to (to,xo) the inverse Lorentz transformation 



7 

-fa 



-fa 
7 



3+1 
-3+1 



t- 



3+1 



(29) 




Fig. 13. A Lorentz transformation acts on the light-cone sim- 
ply as a dilatation. Shown is a discretization invariant under a 
discrete Lorentz transformation, where each point on the plot 
is the Lorentz-transformed of its nearest neighboor towards the 
origin of coordinates. Note this discrete transformation of pa- 
rameter A = j(l + /3) belongs to a subgroup of transformations 
with parameters A 2 ...A n , etc. 



Thus, it appears that we have a simple discretization 
of the line x — t that is invariant under discrete Lorentz 
transformations (dilatations in this line). It is however 
not invariant under translations, as the spacing between 



points Xj 



j > increases in proportion to their 



distance to the origin. Therefore, as is well known, one 
cannot construct a web that is simultaneously translation 
and Lorentz invariant, even under discrete Lorentz trans- 
formations. 

This is different from the Euclidean space case, where 
the translationally-invariant Bravais lattices are also in- 
variant under the discrete rotations of their corresponding 
crystallographic point group. 



Moreover, grids that are invariant under a discrete 
Lorentz transformation, have infinitely many points dense 
at the origin and infinity. 

To construct a discretely Lorentz-invariant lattice of 
the entire 1+1 dimensional space, we just have to write- 
down a discretization in which every point is the Lorentz 
transformed of another (and every Lorentz-transformed 
point belongs to the lattice). Since Lorentz transforma- 
tions leave the metric invariant, the hyperbolae k 2 — (fc ) 2 — 
(k 2 ) = m 2 are invariant, that is, it is sufficient to construct 
a discretization of the hyperbola of (mass) parameter m. 
For this, all one needs to do is to choose a point yo on 
the hyperbola and the discrete Lorentz transformation of 
parameter 7, and apply the Lorentz transformation and 
its inverse to generate the sequence of points 



Vn = A{i n y 



y-n = A(i n y 



Once this has been achieved, it is sufficient to pick up a 
set of hyperbolae to cover the entire space at wish and 
this completes the discretization. 

To construct it, we start with a j/o-lattice in wich there 
are points over both coordinate (x, i)axes. Then we obtain 
two families of hyperbolae with equations 



(30) 



where = ia is the distance between their vertices (which 
will coincide with the chosen yo) an d the coordinate origin. 

Note that the fact that a point and its Lorentz-transformed 
belong to Hi or H2 depends on whether the point lies 
within the forward or backward light-cones (hyperbolae 
of type H2) or not (hyperbolae of type H\) Note also that 
the successive Lorentz-transformed images of a point un- 
der a discrete boost do not fill the hyperbola, they define 
a discrete lattice over it. All that remains is to find this 
subset of points. 

A geometric way to map the discretization of one hy- 
perbola to all others is to simply take straight lines through 
every point on the lattice covering the reference hyperbola 
and the origin. 

This family of straight lines follow either of the rules 
(depending on the sign of the discrete Lorentz transfor- 
mation employed in the construction): 



(ct) = (tanh(i<p))x 
, \ 1 

(ct) = : .. ■ X 

tanh(i^) 



(31) 



where <p is the discrete hyperbolic angle labeling the boost, 
and i the number of direct (positive i) or inverse (negative 
i) Lorentz transformations of the initial vertex. 

Combining eq. (|30|) , and eq. (|3Tj) , we finally construct a 
Lorentz invariant lattice in (1+1) dimmension of a chosen 
tp parameter at the intersections of the two families. The 
situation is plotted in figure [Ml 
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Fig. 14. The intersection points of a family of hyperbolae 
and a family of straight lines define a Lorentz-invariant lattice 
of 1+1 dimensional Minkowski space. The number of points in 
the lattice is infinite for any discrete set of boosts of parameter 
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Now, it's easy to observe that the hyperbolae contain 
all the images under direct and inverse Lorentz transfor- 
mations of the vertex while the family of straight lines 
connect the points that correspond to the A % image of each 
vertex across a family of hyperbolae. 
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